#!/usr/bin/perl
# by lhtk : lhtk80t7@gmail.com
#
# 查找序列中 N 的分布
#
# -------------------------------------
#
# 2019.1.13 改写为通用小工具。
# 
use strict;
use warnings;
use Getopt::Long;

my $fasta;
my $help;
my $f_out;

GetOptions(
    'fasta|i=s' => \$fasta,
    'output|o=s' => \$f_out,
    'help|h' => \$help,
) or die "Error in command line arguments\n";

usage() if ($help);
unless (defined $fasta and -T $fasta) {
    print STDERR "请指定正确的 fasta 文件！\n";
    usage();
}
unless (defined $f_out) {
    print STDERR "请指定输出文件！\n";
    usage();
}
# -----------------------------------------
# 分析
find_Ns($fasta, $f_out);

# ===============================================================
# subroutines
# ===============================================================

sub find_Ns {
    #my $fa = './v1.0/IWGSC_RefSeq_V1_chr1A.fsa';
    my $fa = shift;
    my $out = shift;
    die "$out : exists!" if -e $out;

    my $p = 0; # pos
    my $s;

    open my $fh, '<', $fa or die $!;
    open my $fh_o, '>', $out or die $!;

    my $num = 0;
    while (<$fh>) {
        if (/^>/) {
            $num++;
            if ($num > 1) {
                die "该函数一次只分析一条序列中的 N\n";
            }
            next;
        }
        next if /^\s*$/;

        s/\s//g;
        # 注意一行全是 N，但下一行开始不是 N 了
        unless (/N/i) {
            if (defined $s) {
                print $fh_o $p, "\n";
                $s = undef;
            }
            $p += length;
            next;
        }

        for my $nc (split //) {
            $p++;
            if ($nc !~ /N/i) {
                if (defined $s) {
                    my $e = $p - 1;
                    print $fh_o $e . "\n";
                    $s = undef;
                }
            } else {
                if (not defined $s) {
                    $s = $p;
                    print $fh_o "$s\t";
                }
            }
        }
    }
    print $fh_o "$p\n" if defined $s;
    close $fh;
    close $fh_o;
}

sub usage {
    print <<USAGE;

查找序列中 N 的分布。
注意一次只能分析一条序列。

用法：
    $0 -i input-fasta-file -o output-file

示例：
    $0 -i IWGSC_RefSeq_V1_chr1A.fsa -o chr1A.tab

选项：

    --fasta, -i
        要分析的 FASTA 文件，必须是 DNA 序列，且只含一条序列。

    --output, -o
        输出文件，不可为已存在的文件。
        输出格式为 TAB 格式，每行两列，
        第一列为起始，第二列为终止。

    --help, -h
        打印此帮助。
        可选。

USAGE
    exit 0;
}